doseresNMA antidep/net.structure.R

In htx-r/doseresNMA:

#!!!!!!! I used 200/100 as a doses length but I might make it as an argument
#!!!!!!! I entered 2 the lines for data_no_placebo inside the for-loop but I just test to class=T 
#!!!!!!! I might code the class part to switch drug to class and drug_index to classF instead of repeat the same code 
# Agreed with Georgia: knots <- quantile(min,max)

net.structure <- function(data=data,
                          metareg=F,class.effect=F,
                          ref.lab='placebo',refclass.lab='placebo',
                          cov.pred=0,
                          knot_probs=c(0.25,0.50,0.75)
                          ){
  # 1. data: data.frame() that have the following
    # studyid
    # r: number of events
    # n: sample size
    # dose: the dose level
    # drug: drug names
    # cov (if metareg is TRUE)
    # class (if class.effect is TRUE)
  # 2. metareg: T/F to preform dose-effect NMR
  # 3. class: T/F to run dose-effect class model
  # 4. ref.lab: the referenced drug, default is placebo
  # 5. refclass.lab: the referenced class, default is placebo

  # load
  source('net.functions.R') # it has 3 functions: dose_to_rcs, col_to_mat and direct.comp.index
  require('dplyr')
  require('rms')
  #** initial arguments
  ns <- length(unique(data$studyid))    # number of studies
  ndrugs <- length(unique(data$drug))   # number of drugs
  na <- as.numeric(table(data$studyid)) # number of arms per study
  max.na <- max(na)                     # maximum number of arms

  data$studyID <- as.numeric(as.factor(data$studyid)) # add numeric odered study id
  study_id <- unique(data$studyID)      # a vector of unique studyID

  data$drug_index <- as.numeric(factor(data$drug)) # drug numeric index
  data$drug <- tolower(data$drug) # translate drug names to lower case

   # the numeric index of the reference
  if(is.null(ref.lab)&is.null(refclass.lab)){
    ref.lab<- refclass.lab <-'placebo'
  }else{
    ref.lab<-ref.lab
    refclass.lab <- refclass.lab
  }
  ref.index <- unique(data$drug_index[data$drug==ref.lab])
  refclass.index <- unique(data$classF[data$class==refclass.lab])

  #** doses to RCS transformation

  if(class.effect){ # for hramonized doses
    # dataset without placebo to obtain RCS transformation for each non-placebo drug
    data_no_placebo <- data[data$class!=ref.lab,]
    data_no_placebo$class <- factor(data_no_placebo$class)
    
    nclass <- length(unique(data$class))
    
    rcsdose_drug <- data_no_placebo %>%
      group_by(class) %>%
      group_map(~dose_to_rcs(.x$dose)) # the list order as the order of levels(data_no_placebo)
    rcsdose <- do.call(rbind,rcsdose_drug)
    rownames(rcsdose) <- rep(levels(data_no_placebo$class),
                             times=as.numeric(table(data_no_placebo$class)))
    # ADD rcsdose: dose2
    data$dose2 <- data$dose
    for (i in 1:(nclass-1)){
      data$dose2[data$class==levels(data_no_placebo$class)[i]] <- rcsdose[rownames(rcsdose)==levels(data_no_placebo$class)[i],2]
    }
    
    # dataset without placebo to obtain RCS transformation for each non-placebo drug
    data_no_placebo <- data[data$drug!=ref.lab,]
    data_no_placebo$drug <- factor(data_no_placebo$drug)
    
    # FIND RCS for absolute predictions
    #max.dose <- round(max(data$dose,na.rm = T))
    nd.new <-  100 # predict on 200 doses of each drug
    new.dose <- f.new.dose <- matrix(NA,ndrugs,100) 
    for (i in unique(data_no_placebo$classF)){ # loop through non-placebo drugs
      max.dose <- round(max(data[data$classF==i,]$dose,na.rm = T))
      min.dose <- round(min(data[data$classF==i,]$dose,na.rm = T))
      new.dose[i,] <-  seq(0,max.dose,length.out = 100)
      knots <- quantile(min.dose:max.dose,probs = knot_probs)
      f.new.dose[i,] <-  rcspline.eval(new.dose[i,],knots = knots,inclx = TRUE)[,2]
    }

  }else{ # for unhramonized doses
    # FIND rcs
    rcsdose_drug <- data_no_placebo %>%
                     group_by(drug) %>%
                      group_map(~dose_to_rcs(.x$dose)) # the list order as the order of levels(data_no_placebo)
      rcsdose <- do.call(rbind,rcsdose_drug)
      rownames(rcsdose) <- rep(levels(data_no_placebo$drug),
                               times=as.numeric(table(data_no_placebo$drug)))
    # ADD rcsdose: dose2
    data$dose2 <- data$dose
    for (i in 1:(ndrugs-1)){
      data$dose2[data$drug==levels(data_no_placebo$drug)[i]] <- rcsdose[rownames(rcsdose)==levels(data_no_placebo$drug)[i],2]
    }
    # FIND RCS for absolute predictions
    #max.dose <- round(max(data$dose,na.rm = T))
    nd.new <-  200 # predict on 200 doses of each drug
    new.dose <- f.new.dose <- matrix(NA,ndrugs,200) 
    for (i in unique(data_no_placebo$drug_index)){ # loop through non-placebo drugs
      max.dose <- round(max(data[data$drug_index==i,]$dose,na.rm = T))
      min.dose <- round(min(data[data$drug_index==i,]$dose,na.rm = T))
      new.dose[i,] <-  seq(0,max.dose,length.out = 200)
      knots <- quantile(min.dose:max.dose,probs = knot_probs)
      f.new.dose[i,] <-  rcspline.eval(new.dose[i,],knots = knots,inclx = TRUE)[,2]
    }

  }
  # ORDER the doses on each study: start from zero
  data_withRCS <-data%>%
                  group_by(studyID)%>%
                   arrange(dose,.by_group=TRUE)

  #** Column to  matrix with ns in rows and na in columns(max.na is specified, the additionals are given NA)
  rmat <- col_to_mat(data_withRCS,data_withRCS$r)

  nmat <- col_to_mat(data_withRCS,data_withRCS$n)

  dosemat <- col_to_mat(data_withRCS,data_withRCS$dose)

  dose2mat <- col_to_mat(data_withRCS,data_withRCS$dose2)

  tmat <- col_to_mat(data_withRCS,data_withRCS$drug_index)

  #** For inconsistency model: indices of direct comparisons
  comp <- direct.comp.index(data_withRCS)
  t1 <- comp[,'t1']
  t2 <- comp[,'t2']
  ncomp <- nrow(comp)

  #** RETURN: list to jags models
  jagsdata <- list(ns=ns, na=na,ndrugs=ndrugs,
                                 r=rmat, n=nmat,dose1=dosemat,t=tmat,dose2=dose2mat,
                                 ref.index=ref.index,
                                 rr=rmat,nn=nmat,new.dose=new.dose,f.new.dose=f.new.dose,nd.new=nd.new,
                                 t1=t1,t2=t2,ncomp=ncomp
                                 )
  #** Add to jagsdata for metareg or class effect parts 
  if(metareg){
    if(class.effect){
      classmat <-col_to_mat(data_withRCS,data_withRCS$classF)
      nc <- length(unique(classF))
      covmat <-col_to_mat(data_withRCS,data_withRCS$cov)
      add <- list(class.effect=classmat, 
                  refclass.index=refclass.index, 
                  nc=nc,
                  cov=colMeans(t(covmat),na.rm=T),
                  cov.pred=cov.pred
                  )
      jagsdata <- c(jagsdata,add)
    }else{
      covmat <-col_to_mat(data_withRCS,data_withRCS$cov)
      add <- list(cov=colMeans(t(covmat),na.rm=T),
                  cov.pred=cov.pred
                  )
      jagsdata <- c(jagsdata,add)
      }
  }else{
    if(class.effect){
      classmat <-col_to_mat(data_withRCS,data_withRCS$classF)
      nc <- length(unique(data_withRCS$classF))
      add <- list(class=classmat, refclass.index=refclass.index, nc=nc)
      jagsdata <- c(jagsdata,add)
    }else{
      jagsdata
    }
  }

  return(jagsdata)
}